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Methods of modeling cellular regulatory networks as diverse as differential equations and Boolean 
networks co-exist, however, without any closer correspondence to each other. With the example 
system of the fission yeast cell cycle control network, we here set the two approaches in relation to 
each other. We find that the Boolean network can be formulated as a specific coarse-grained limit of 
the more detailed differential network model for this system. This lays the mathematical foundation 
on which Boolean networks can be applied to biological regulatory networks in a controlled way. 
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I. INTRODUCTION 

The task of molecular cell biology is to comprehend the control of cellular processes of living cells encoded in 
the genome of the cell. These cellular processes are guided by sophisticated networks of interactions between 
the macromolecules of the cell as proteins, nucleic acids and polysaccharides. Their complexes and structures 
define the unique properties that enable them to perform the functions of the cell as, for example, catalysis of 
chemical transformations, production of movement, and heredity. The complexity of these processes demands 
not only advanced experimental techniques, but also adequate mathematical and computational models for 
understanding them (Gunsalus et al., 2005; Riel, 2006). 

Today, there are different methods for modeling the complex networks of biochemical interactions, ranging 
from master equations based on first principles (Monte-Carlo method) (Gillespie, 1976; Gillespie, 1977), 
ordinary differential equations (ODE) (Aguda, 2006; Chen et al., 2000; Novak and Tyson, 1993; Novak and 
Tyson, 1997; Novak et al., 2001; Novak and Tyson, 2004; Sveiczer et al., 2000; Tyson et al., 2001; Tyson et 
al., 2003), stochastic differential equations (Fokker-Plank equations), all the way to Boolean networks (Albert 
and Othmer, 2003; Bornholdt, 2005; Sanchez and Thieffry, 2001; Thomas et al., 1995). 

Among these methods, a most popular approach to modeling biochemical networks is via differential equa- 
tions, based on the known chemical kinetics and successfully applied to describing numerous processes in living 
organisms (Chen et al., 2000; Novak and Tyson, 1993; Novak and Tyson, 1997; Novak et al., 2001; Novak and 
Tyson, 2004; Tyson et al., 2002). To build an ODE model, one starts with a schematic diagram representing 
the known interactions between components. Then this diagram is converted into a set of differential and 
algebraic equations. The full ODE model then consists of this set of rate equations, plus a set of parameter 
values as well as a set of initial conditions. The solutions of the ODEs give the time-dependence of each 
component of the system. In practice these solutions depend on rather detailed knowledge about all reactions 
and kinetic parameters. 

In studies where prediction of exact reaction times is not an issue, simpler models than ODE and less 
parameters may be necessary for predicting the course of events in a regulatory network. For example, 
relevant features of cell commitment, cell cycle progression, and cell differentiation are already described 
in terms of a sequence of regulatory events (Albert and Othmer, 2003; Braunewell and Bornholdt, 2006; 
Davidich and Bornholdt, 2008; Li et al., 2004; Sanchez et al., 2001). In such cases, the much simpler modeling 
framework of Boolean networks may be a suitable method (Kaufmann, 1969; Thomas, 1973). Constructing 
a Boolean network model starts from a wiring diagram of interactions between biochemical elements as well, 
but no kinetic details are needed. Interactions are classified into just two classes, activation or inhibition, as 
well as the concentration levels being reduced to just an ON or OFF state. 

Despite their extreme simplicity, such Boolean models are able to reproduce regulatory sequences, for 
example, models of genetic networks of A. thaliana, (Espinosa-Soto et al., 2004; Mendoza et al., 1999; Thum 
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et al., 2003), the cell-cycle network of S. cerevisiae (Li et al., 2004), the mammalian cell-cycle (Faure et al., 
2006), and the segment polarity gene network in D. melanogaster (Albert and Othmer, 2003; Sanchez et 
al., 2001). These examples show that the Boolean network approach provides reliable results for different 
organisms. 

The two diverse methods just summarized (ODE and Boolean networks) are both based on the same 
"wiring" diagram of interactions between the components, however, use much different amounts of information 
about these interactions. This poses the interesting question how these two methods are related to each 
other. A first correspondence between Boolean network and ordinary differential equations has been drawn 
by Glass and collaborators (Glass and Kauffman, 1973; Glass and Hill, 1998) who explored the relationship 
between a class of non-linear equations representing biochemical control networks and homologous switching 
networks. They argued that such a correspondence can be achieved with the following requirements: 1) Rates 
of reactions are described by monotonic sigmoidal functions having distinct upper and lower asymptotes. 2) 
The parameters must be defined to match the upper or lower asymptote. 3) The target control function must 
correspond to the maximal or basal rate of biochemical processes. They subsequently demonstrated that there 
is a large variety of such functions as, e.g., the Heaviside function, the error function, or the Hill function 
defined for positive arguments. 

This leads to a mapping between asymptotical solutions of the ODE system and the Boolean system, while 
omitting the exact way of transitions between dynamical states. 

In this paper, we further explore the correspondence between ODE and Boolean network models considering 
a specific biological system and demonstrate how a working Boolean model can be derived in terms of a 
mathematically well defined coarse-grained limit of an underlying ODE model. As our working example 
we choose the fission yeast cell-cycle control network (Schizosaccharomyces Pombe). The division cell-cycle 
consists of four phases G1-S-G2-M, during which DNA is replicated and the cell divides itself into two cells. 
The main role is played by " Cyclin-Dependcnt-Kinascs" (CDKs) and cyclins that bind to CDKs to form 
complexes. CDKs, while being present at all times, can only be active in complexes with cyclins. Cyclins 
are synthesized or degraded depending on other regulatory activities. Another important participant of the 
process is an enzyme complex called the " Anaphase-Promoting Complex" (APC), which targets cyclins for 
degradation. In summary, "to understand the molecular control of cell reproduction is to understand the 
regulation of CDK and APC activities" (Tyson et al., 2002). These processes, while being complicated, have 
been well studied for fission yeast S. Pombe and successful ODE models (Novak et al., 2001) exist. We choose 
the most widespread version of the model (Novak et al., 2001), as often cited in textbooks, which will serve 
us as a starting point for our study. 

The paper is organized as follows. In Section 2 we show the passage from the ODE system of algebraic 
differential equations for the fission yeast cell cycle, to the limit of the corresponding Boolean model which 
we construct. Here also the difficulties that one can meet when works with Boolean approaches are discussed. 
Section 3 explores the dynamics of the derived Boolean model of the fission yeast cell-cycle. 

Finally, in the discussion section the properties of the obtained system are recapitulated and the Boolean 
and ODE approaches are compared. 

II. BOOLEAN VARIABLES AS STATIONARY STATES 

The passage from a differential equations model to a Boolean model requires the mapping of continuous 
solutions (Novak et al, 2001) into the ON/OFF states of a Boolean network's nodes. In order to achieve this, 
the time evolution of a function, determined by the rate functions and kinetic constants, has to be replaced 
with a discrete mapping of the node set into itself. Moreover, the rules of this self-mapping have to be governed 
by logical functions, connecting the binary states of interacting nodes. The dynamics of the resulting Boolean 
network is an ordered sequence of states of the network nodes, instead of the continuous time output of the 
ODE model. Having this in mind, let us find the conditions, which allow to perform the transition from a 
differential equations model to a Boolean system. In the following, we will first describe the passage of the 
continuous variables to discrete states and, in a second step, construct the logical functions representing the 
dynamics. 

A. Stationary states of ODE system 

The model of the fission yeast cell-cycle (Novak, 2001) is based on the antagonist interaction of CDK- 
cyclin complex with APC (via proteins Ste9 and Slpl) and CDK inhibitor Ruml. The CDK-cyclin complex 
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(Cdc2/Cdcl3) is represented by two variables - preMPF and MPF ("Maturation Promoting Factor"). 
Furthermore, helper molecules (Start Kinase SK, transcription factor TF, kinase k wee , phosphatase ^25, and 
time-delay enzyme IE) participate in the process. Regulatory interaction between these macromolccules are 
described by two types of equations - differential equations for Ste9, Slpl, IEP, M, SK, Rural, preMPF 
and algebraic equations for k wee , k 2 5, TF and MPF variables. Among the first ones, some are of Michaclis- 
Mcnten type (Ste9, Slpl, IEP) and some are exponential growth (M, SK) equations. The mass variable M 
plays a special role as it parametrizes the time evolution in the system. The ODE model (Novak, 2001) uses 
arbitrary units for concentrations in all equations, since there are few data of actual protein concentrations. 
The kinetic constants deretmine the right timing of the processes. Solutions of the system show that the 
concentrations of the major proteins in general rise or decrease steeply. 

To make the transition to a Boolean system, we first need to rescale the differential equations such that their 
solutions assume values between (inactive) and 1 (maximum activity). This is a first step towards mapping 
these variables onto Boolean OFF / ON variables with values and 1 . The rescaling does not change the form 
of equations, it only affects the values of kinetic constants. To do this, let us divide all functions by their 
respective maximum value. For example, for Slpl we introduce the new rescaled function Slpl\=SlplAmpl, 
where Ampl = 2.1 is the amplitude of the original solution. Rescaling all variables except M we obtain: 
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where the Goldbeter-Koshland (GK) (Goldbeter and Koshland, 1981; Novak et al., 2001) function has the 
following general form: 



G(a, b, c, d) = 



2ad 



b- a + bc + ad+ y/ (b - a + be + ad) 2 - Aad(b - a) 



(16) 



Square brackets denote the concentrations of their elements. The subscript 1 marks the rescaled variables 
with maximum 1. The new values of parameters are shown in Table 1. 

Next we map the continuous solution of the ODE model (Novak et al, 2001) into the discrete states of a 
Boolean network's nodes. Since in a Boolean model there is no continuous time, but rather a sequence of 
switching events between two stationary states of the nodes, one needs to reduce (whcrccvcr possible) the 
initial dynamics of the system to a sequence of the evolution of stationary states. For this, one can use the 
results of a bifurcation analysis (Novak et al., 2001) of the transitions during the cell cycle. Thus, there are 
some variables (Ste9, Slpl, IEP) that are described by Goldbeter-Koshland (GK) functions (Novak et al., 
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2001) in the stationary state: 

[Ste9i] = G(k' 3 + k^[Slpl 1 ],k' 4 [SK 1 ] + k 4 [MPF 1 },J 3 ,J i ) (17) 

[IEPi] = l/k' 9 G{k 9 [MPFi],kio,J9,Jio) (18) 

[Slpli] = [Slpl T i]G(kr[IEPi],k8,J7/[SlplTi],J8/[Sl P l T i]) (19) 

The characteristic properties of the GK function imply that its variable mainly resides in two limiting 
states: High and Low. The transition time between them is short (as the variables J are small). Therefore, 
we approximate them as Boolean (binary) variables 



It is easy to see that Slplxi determines only the amplitude and the smoothness of the transition in (191, 
therefore, we neglect Slplxi as a first step towards a Boolean model and write 

[Slpli] = G(krIEPi,k a , J, J). (20) 

The states of the remaining variables can be described by Golbeter-Koshland functions, as well. For SK, 
while the equation for this variable is exponential, one can evaluate the right-hand part through a GK function 
with the stationary solution (Novak et al., 2001): 

[SKi] = {kis/ku)[TF] = {k ls /k 14 )G(k 15 M 1 ,k' 16 + k'{ 6 [MPF 1 ],J 15 ,J 16 ). (21) 



Furthermore, there are three algebraic equations for TF (10 1, k wee (11), and ^25 (12 1 that contain Golbeter- 
Koshland functions. Here, again, the two limiting states of the GK function will be related to the binary 
ON/OFF version of the corresponding variables in the Boolean limit. 

Finally, we will simplify the functional behavior for the CcIcISt, preMPF, and Rumlx, as well. Again we 
want to neglect the exact path of their transitions, keeping the limiting stationary states, eventually enabling 
us to take the limit of Boolean functions as a simplified description of the dynamics. Equations with the 
following requirements will allow us to take this limit: 

a) In a small neighborhood of the switching point, the functional behavior can be approximated by an 
exponential rise. 

b) On the larger interval, it has a stationary solution with the steep transition between two limiting 
stationary states. 

c) This function converges to the Heaviside function in the limit of steep transition. 

In the following we will see that the Michaelis-Menten dynamics fulfils these requirements, resulting in exact 
conformity of initial and final states and permitting a well controlled passage to a Boolean function. Let us 
start from the Michaelis-Menten equation 

dX , l-X , X 

¥ = il j^n-fe (22) 



and first check the condition a). Expanding (22) in the neighborhood of switching points where X « J x 2 << 
1 and keeping leading order terms yields 

This is a common equatio n of exponential growth/decrease. This allow us to take equations of exponential 



growth as an expansion of (23 1 in the neighborhood of switching ON/OFF points. For CdclSx, for example, 
the equation 

d \ Cdc ^ T ^ = klM - (k' 2 + k%[Ste9i] + fc2"[SZpli])[C*dcl3 T i] (24) 
is the expansion of the equation 

dJ£ ^ - - « + + t:i *« m - | J , [ "Xr (25) 

For illustration, in Fig. 1, this function is compared to the initial CdclSxi- 



Condition b) is satisfied as well, since the stationary states of (22 1 are described by a Goldbeter-Koshland 
function. Validity of condition c) is shown in the next section. 

Thus, for the above equations we have a system of GK-functions which are responsible for the transitions 
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between stationary states: 



[Cdcl3ri] = G(kiM,k%[Ste9i] + k2[Slpli],J,J) (26) 
\preMPFi) = G(Mw e [Cdcl3 T i], [k weei ] + fco'[fc 2Sl ] + k' 2 + k£[Ste9i] + k! 2 "[Slpl x ]),J,J) (27) 
[Rumlri] = G(fcn, fci2 + k' 12 [SKi] + k" 2 [MPFi], J, J). (28) 

One can see in Fig. 1 that the new and initial functions start to grow and start to decrease at the same times, 
respectively. Note that the obtained substituted equations ( 26 28 1 play only a helper role and cannot be 
directly applied to the initial system of differential equations. The complete correspondence of the obtained 
system to the initial one is achieved by the limit transition shown in sectio n |IIB| 

Summarizing all information above we obtain the system of equations ( 10|12 1, ( 17|l8 ), ( 20pT l, ( 26]|28 1, 
describing the stationary states of the corresponding variables. The next step is to perform a transition 
from the continuous functions to discrete functions. In functions with continuous time, the stepwise change 
corresponds to the Heaviside function. 



B. Passage to Boolean variables 



We now show that there is an exact passage from the function of Goldbeter-Koshland (Goldbeter and 



Koshland, 1981) to the indicator function (Heaviside step function). Let us remark that in ( 16 1 the parameters 
a and b are functional variables, whereas c and d (in ( To] 12 1 , ( T7fl8 l, ( 20pT ), ( 26p8 1 they are denoted as 
variables J) are usually fixed and small in all equations. The range of values of c, d varies from 0.001 to 0.01 
and only one time for Slpl it takes on the value 0.3. Very small parameters c and d mean that the enzyme- 
substrate complex is tightly bound and hardly dissociates. Thereby in (Novak et al., 2001) an assumption is 
made that the enzyme-substrate complexes involved are very stable (Aguda, 2006). For this reason, let us 
consider the behavior of the corresponding GK-functions in the limiting case c — > 0, d — > while a and b take 
finite values. 

First note that, both, numerator and denominator depend on d. Moreover, a is a numerator's factor. At 
the same time, c appears in a sum with the finite terms in the denominator, only (except at the point b = a). 
Therefore, we can assume, without loss of generality, that c = and consider below 



G(a,b,0,d) 



2ad 



b — a + ad+ *J (b — a + ad) 2 — Aad{b — a) 



Using the Taylor expansion of the square root in the denominator and neglecting higher powers of d, we 
obtain for all points a =/= b 

G(a,b,0,d) = - , ■ 2 , Qrf -s— . (29) 



b — a + ad + \b — al — ad-^— % 

I I |d— a\ 



There are two possible cases: 



a) if b — a > then \b — a\ — b — a and ( 29 ) takes a form 



G(a, b,0,d) 



2nd 



b) if b- a < then \b ■ 



(b-a) and (29) is simply G(a,b,0,d)=l. 



This implies that there is a passage to the limit from Goldbeter-Koshland function G(a, b, c, d) to the Heaviside 
function 8(a — b): 

' 0, a < b, 



lim 



?(a-6) 



1, a > b. 



Thus, in this limit, the variable a plays the role of an activator input and the variable b is an inhibitor input. 
The output is active/inactive (its Boolean value is equal to one/zero) if the total value of activator inputs is 
larger/smaller than the total value of inhibitor inputs. Thus, the Goldbeter-Koshland-function converges to 
the Heaviside function in the limit of steep transitions. 



G 



C. Logical Boolean functions 

Let us now rewrite the system of equations in the limit of the parameters J — + ( J3, J5, J7, Jg, JiO, Jig, 

Jawee, Jiwee, Ji25, Ja25 and the corresponding parameters J in equations (20), (26-28) equations) as 

[SteQi] = 9{k' 3 + k'slSlpli] - k'^SKi] - k^MPFi]) (30) 

[IE Pi] = 9(k 9 [MPFi]-ki ) (31) 

[SKi] = (ki 3 /ki4)9{kisMi-k'i 6 -k'UMPFi]) (32) 

[TFi] = e(ki 5 M-k'i 6 -ki 6 \MPFi]) (33) 

[^ii?eei] = ^ruee {^wee k wee )9{\^ awee Viwee [MPF± ] (^4) 

[fcasj = k' 25 + (k'l B -k' 25 )e(y a 25[MPFi]-V i25 ) (35) 

[Cdcl3Ti] = 9(kiM -k 2 lSte9i]- k' 2 [Slpli]) (36) 

[preMPFi] = 8(k wee [Cdcl3 T i] - [k weei ] ~ [k 25l ] - k' 2 - k' 2 \Ste9i] - k' 2 "[Slpli]) (37) 

[RumlTi] = e(ku-ki 2 - k'i 2 [SKi]-ki 2 [MPFi] (38) 

[Slpli] = e(k 7 [IEPi]-k s ). (39) 



Let us add two simplifications to the equations (32 1 and ( 34p5 1. In equation (32 1, the coefficient fcis/fcu = 1 



and thus can be neglected. In equation (34 1, k weei can have two possible values: 0.115 and 1. The first one 
can be reduced to since it does not change the behavior of the system and analogous for fc 2 5 

[^weei] 9(V a wee> Viwee [hi P F\] , Jawee, Jiwee) (40) 

[k 25l ] = 9([MPFi]-V i25 ) (41) 
[SKi] = 9(k w M-k'i 6 -kUMPFi]). (42) 

As a result we have a system of ten equations, where all variables, except MPF and M take values or 1. 
MPF and M cannot be described in this formalism. MPF is represented by an algebraic equation which 
cannot be reduced to the GK- function. Taking a closer look, the solution of MPF does not reach a simple 
stationary state, instead there are three typical states of MPF in the system (preMPF, Ruml, Cdcl3) - 
OFF, intermediate and high activation. 

1) If Ccfcl3 = then MPF = 0, independently of the states of preMPF and Ruml. 

2) If Cdcl3 = 1 and preMPF = 1 then MPF — 0.14, independently of the state of Ruml. This corresponds 
to an intermediate level, where preMPF prevents high excitation. 

3) If Cdcl3 = 1, preMPF = 0, and Ruml = 1, then MPF = 1, with its value slightly decreasing to 
MPF = 0.93 if Ruml = 0. This corresponds to a high level of activation, when MPF is activated by Cdcl3 
and this activation is not reduced by preMPF. 

Let us reformulate these rules in the following. Assume there are two variables - MPF and MPF2. The 
first one is activated by Cdcl3. For activation of the second variable MPF2, one assumes that MPF has 
to be present as a low-level and preMPF should be inactive. Thereby, one needs to rewrite the system of 
equations (30 1, ( 36|37 ) taking into account which level of excitation of MPF is crucial for each particular 
variable. 

(43) 
(44) 
(45) 
(46) 
(47) 
(48) 
(49) 
(50) 
(51) 



[MPFi] 


= 6(Cdcl3 T i) 


[MPF2] 


= 0(MPFi - [preMPF!]) 


[IE Pi] 


= d([k 9 [MPFi] - kio) 


[TFi] 


= 9(k 15 M -k'i 6 -k , ; fi [MPF2 1 ]) 


[RuraXn] 


= 9(kn - fc u - k[ 2 [SKi] - k'{ 2 [MPFi]) 


\kwee± ] 


= 0(V awee -V ituee [MPFi]) 


mi 


= 0{[MPFi]-V i25 ) 


[SKi] 


= 0{[TF]) 


[5*691] 


= 0(k' 3 + k'^SlpU] - ki[SKi] - k 4 ([MPFi])) 



Second, as in the model based on differential equations the cell mass M takes a special role in the present 
model. The solution (Novak et al., 2001) treats it during a cell growth as an independent variable, which is 
described by an exponential growth function Q. Thus, the variable M corresponds to a time in this system, 
which drives the evolution between stationary states (Novak et al, 2001). In the system, M directly influences 
Cdcl3 and TF. As soon as M reaches a threshold value, it activates Cdcl3 and induces the sequence of 
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consecutive transitions between stationary states. For TF it plays a role of constantly positive input, TF is 
always active unless MPF has a high activity. 

As a criterium for the end of the cycle, Novak et al. (Novak et al., 2001) determine when the cell divides by 
monitoring the values of the other variables. When these chosen variables have certain values that indicate 
the end of the cell cycle, the current value of M is divided by two manually, as at the end of mitosis the cell 
divides into two daughter cells of approximately equal masses. Subsequently, M continues its exponential 
growth, again. 

Following this strategy, one needs to distinguish M between two principal different values - M and 2M in 
the Boolean model. Here M works at the beginning of the cell cycle as a trigger of switching events, whereas 
2M play a role of an indicator for the end of the cell cycle. Correspondingly, M becomes 2M at the end of 
mitosis, when Slpl, Ste9 and IEP all have high concentrations. 

Thus, one can add the following Boolean rule: 

M = 6(2M- [Ste9i][Slpli][IEPi]) (52) 
2M = 6{M[Ste9i][Slpli][IEP x ] -2M). (53) 



Thus, we have a system of equations ( 43||53 1, where each variable can take values or 1, only. It is easy to 
simplify this system, reducing the kineti c co emcients to or 1 and adding thresholds. Consider, for example, 
Cg?c13t- In Table 2, based on equation (36), we show all possible cases for Cdcl3. In a more compact form, 
where the kinetic constants are reduced fol, these rules become 



[CdclSTi] = 6(M - [Ste9i] - [Slpli]). 



(54) 



Repeating the same procedure for all variables, we obtain the system of equations: 



[preMPFi] = 6(k wee + [Cdcl3 T i] - 1 - [fe 6l ] - [Ste9i] - [Slpli])) (55) 

[Slpli] = 6([IEPi]) (56) 

[TFi] = = 6([M] + [2M] - [MPFIx] (57) 

[IEP!] = 9([[MPF 2 ]) (58) 

[RumlTi] = 6(0.5 — [SKi] — [MPFi]) (59) 

[k weei ] = 61(0.5 - [MPF!]) (60) 

[fe 6l ] = 0{[MPFi] - 0.5) (61) 

[SK!] = 9([TFi]) (62) 

[Ste9i] = 0([Slpli] - [SK!] - [MPF!]) (63) 

M = <9(2M + 3 - [Ste9i] - [Slpli] - [IEP!]) (64) 

2M = 0(M + [Ste9i] + [Slpli] + [IE Pi] - 3 - 2M) (65) 



III. BOOLEAN MODEL 



We now have a system of algebraic equations (53631 
stationary states, plus Boolean equations ( 43|44 ),( 64|65 
no information about continuous time is present any more, 
discrete dynamical sequence, we iteratively solve the system. 
(Novak et al., 2001): 2M = 1, Slpli = 1, IEPi=l, Stt&x = 1, 



which describe the switch-like transitions between 
for M and MPF. Note that in this discrete system, 
except the sequence of events. To obtain this 
We start from the known initial conditions 
kweei — 1, with all other variables being 0. 
Following the terminology of Boolean models, each variable is represented by one node. The network of nodes 
is shown in Fig. 2. Each node i has only two states, Si(n) = 1 (active) and Si(n) = (inactive). The index n 
is the number of iterations. The iterative solution of the system has the following general form: 



Si(n+l) = 9[Y,T ik S k (t) 



(66) 
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with the Heaviside function 

*(*)={ J; *lo° ^ 

All nodes are updated synchronously, which corresponds to the iteration of the full dynamical system. The 
interaction matrix and the state vector Qi determine the transition rule between states. 

The specific values of the interactions are determined as follows. All interactions between a pair 
of nodes are defined by an interaction strength given as the elements and an activation threshold Qi. 
Positive (negative) arguments of the ^-functions have = 1 (T^ = — 1). Tu = 1 for M, this rule is true 
only for M node, since it is described by a growing exponential function. 

The resulting matrix and vector Qi have the following form: 
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IV. RESULTS OF BOOLEAN SIMULATION OF THE FISSION YEAST CELL CYCLE 



First we run the model described in the previous section with the initial conditions (Novak et al., 2001). 
The temporal evolution of the protein states is presented in Table 3. One can see that the iterative solution 

) is the switching between unstable stationary states, which coincide with the 
corresponding evolution in the ODE model. The final state is a stable stationary state of the system. One 
notices that the initial and end states are identical except for the activation of the nodes M and 2M. The 
update of nodes M and 2M, keeping all other nodes in the same states, starts the new cycle. This cycling of 
the model is similar to the realization of cycling in the original model of differential equations. 
Let us briefly summarize our coarse-graining strategy that we followed in this article. 

In Fig. 3, we show consecutive abstractions of the model for the Ste9 and IEP variables as an example. For 
this we first plot the dynamics of the differential equations model, then the evolution between the stationary 
states solutions of the ODE model, and finally the sequence of states obtained from the iterated Boolean 
network model. 

In the next step we run the model starting from each of the 2 15 possible initial states. We find that from 
all initial states 67% flows into one big attractor. This attractor is the same stable stationary state that one 
obtains starting with biological conditions described above. 



of the system (43 44) (54p5 



A. Mutations 



Let us also compare the behavior of mutants. We model two mutations - Wee~ and Wee~Cdc25A described 
in (Novak et al., 2001). In Boolean models one cannot distinguish between reduced activity and no activity. 
This is why we model Wee~ as WeeA in both cases. We run the model starting with wild type initial 
conditions, but with deleted nodes k wee , and in the second case k wee and k 2 5- In both mutations, the number 
of steps is reduced to 12, as compared to 14 in the wild type cell cycle described in the previous section. This 
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suggests that the cell can divide at a smaller size than the wild type, where both mutations are viable. Our 
results are in accordance with the predictions of the earlier differential equation model (Novak et al., 2001). 

B. Comparison with an existing Boolean model for the fission yeast cell cycle 

It is interesting to compare this model with another Boolean model for the fission yeast cell cycle (Davidich 
and Bornholdt, 2008) that was built on known biochemical interactions between proteins, only. Both models 
are quite similar and have the same connections between homologue nodes. Their dynamics matches the 
wild type signaling sequence during the cell cycle. The difference is that in the current model the nodes 
Cdcl3, preMPF, MPF1, MPF2 correspond to only two nodes Cdc2/Cdcl3 and Cdc2/Cdcl3* in (Davidich 
and Bornholdt, 2008). So in the model (Davidich and Bornholdt, 2008), the complex Cdc2/Cdcl3 can have 
two levels of activation - medium and high. The intermediate level corresponds to sole activation of the 
Cdc2/Cdcl3 node, whereas a high level of activation is represented by activation of both, Cdc2/Cdcl3 and 
Cdc2/Cdcl3* .The last one plays the role of a dephosphorylated Cdc2/Cdcl3, that closely corresponds to 
MPF in (Novak et al., 2001). In the current model, the node MPF was separated into two nodes MPFl 
and MPF2 as well, to distinguish different levels of MPF activation. 

The formulas responsible for evolution of proteins are similar in a sense that in both a threshold activation 
rule is used. In the model (Davidich and Bornholdt, 2008), proteins remain active if the corresponding node 
was not switched-off by other incoming inhibiting signals. This rule means that if the protein was activated 
it requires some other signals to change its state. Whereas in the current model one needs to have always a 
positive incoming signal in order to keep the protein in its active form. 

V. DISCUSSION AND CONCLUSION 

In this paper, our aim is to make a connection between two successfully used methods - the ODE and the 
Boolean models for predicting properties of a real biological process. In particular, we show a possible limit 
transition between the ODE and the Boolean model for the fission yeast cell cycle. For this purpose the known 
ODE model (Novak et al., 2001) of fission yeast has been chosen. 

In order to do the transition, first the differential equations are rewritten in a such way that the solutions 
of equations reach 1 in their maximum values. Then the obtained equations are transformed in a limiting 
procedure to Boolean functions. Firstly, Michaelis-Menten equations and equations with switch-like behavior 
can be directly reduced to Boolean functions. Secondly, a set of equations with sigmoidal transfer functions 
can be replaced with Michaelis-Menten equations without changing the sequence of states through which 
the system evolves. Thirdly, there are also some cases that cannot be reduced to the two previous ones. It 
happens when the variable is described by a constantly growing function or a function which has distinctly 
different levels. In this case we propose in the Boolean model to substitute those variables with two labeling 
intermediate and high activity of it. Finally, all continuous solutions of equations are mapped into ON/OFF 
states of Boolean network and the transition between states are described by Boolean functions. 

This Boolean model reproduces the results of the initial ODE model (Novak et al., 2001). In particular, 
starting with initial conditions as in (Novak et al., 2001), the system evolves through the same sequence 
of states. The second evidence of similar behavior of the ODE and the Boolean model is the robustness 
to the initial conditions. The Boolean model has a dominant attractor (67%), attracting most of the tra- 
jectories, starting from different initial conditions. The dominant attractor coincides with initial biological 
conditions of the system. The ability to model mutations in the Boolean approach additionally confirms a 
good correspondence between the ODE and the Boolean model. 

We find that the transition to a Boolean model is possible for differential equations, which have monotonic 
sigmoidal functions with distinct upper and lower asymptotes. In particular, firstly, in our case Michaelis- 
Menten equations are reduced to S-shaped GK-functions which have the necessary asymptotes (Glass and 
Kauffman, 1973; Glass and Hill, 1998). This function works as a switch between the cases when parameters 
are defined as the upper or lower asymptote and the target control function corresponds to the maximal or 
basal rate of biochemical processes. Secondly, here substituting some equations that have monotonic sigmoidal 
functions on the right-hand-side with Michaelis-Menten functions, we also find, that the exact form of the 
sigmoidal function does not strongly influence the behavior of the system. The comparison of the current 
model with a previous Boolean model for fission yeast reveals that they both have a similar set of variables 
(proteins) and similar Boolean functions responsible for update. 
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Our results also confirm the idea that some molecular control networks are so robustly designed that timing 
is not a critical factor (Brauncwcll and Bornholdt, 2006). In our case it is possible to reproduce the main results 
of (Novak et al., 2001) without including time, but reproducing the right sequence of events. It supports the 
idea that the Boolean approach could contain sufficient information. Thereby one needs less information about 
the system, the knowledge about reactions on the level of activation/inhibition is sufficient, which eliminates 
the problem of finding the right kinetic constants. Another advantage is the low computational cost of Boolean 
networks. The problems one meets working with the Boolean approach are that it is sometimes difficult to 
reduce the concentration level of some proteins only to ON/OFF states. Sometimes there are intermediate 
states of concentration which need to be separated from high concentration. In this case two methods are 
possible. One is, as we implemented it here, to divide this variable into two and to perform as two different 
nodes in a system. Doing this, one needs to take into account the differences in influences of this protein when 
it has intermediate and high concentration. Another solution for a such situation could be the introduction of 
two discrete levels of concentration that the protein can have, for example 1 for intermediate and 2 for high 
concentration, which was used in modeling gene regulatory network model for Arabidopsisthaliana flower 
development (Espinosa-Soto et al., 2004). 

We would like to note that the ODE and the Boolean approach are both useful methods. The advantage of 
the ODE approach is that it provides detailed information about the system at any given time in contrast to 
the Boolean method, which reproduces only the right sequence of events. But the costs for this information 
are the following. One needs to have exhaustive information about the reactions, where the most difficult 
part is to find the right kinetic constants. Also it demands more computational costs to find the solution 
of the system. One could say that the ODE approach is appropriate when the system is well studied and it 
is necessary to make a detailed study of all reactions that take place. On the other hand if the task is to 
understand the main principles of some process and one has less information, the Boolean approach is very 
suitable to use. 
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FIG. 1: Numerical simulation of (24 1 (black curve) and (25 1 (red curve). Note that the downward decrease occurs 
practically simultaneously. 



Cdcl3i 


fci = 0.04, k' 2 = 0.03, k% = 1, fc 2 " = 0.21 


preMPF 1 


k' = 1.5,fc£ = 1.17, fc'"0 = 5 


Sie9i 


fc^, = 1, fc£ = 21, J 3 = 0.01, fc 4 = 1.98, k 4 = 50.75 


Slpl T1 


k' 5 = 0.002, fc£ = 0.143, k' 6 = 0.048, J 5 = 0.20689 


Slpli 


k 7 = 0.429, k 8 = 0.119, J 7 = 0.0005, J 8 = 0.0005 


IEP 


fc 9 = 0.16, J 9 = 0.01, fcio = 0.01, Jio = 0.011 fc 9 = 0.91 


Rural 


feu = 0.698, fci2 = 0.01, k' 12 = 0.99, k'{ 2 = 4.35 


SK 


fcl 3 = 0.1, fcu = 0.1 


M 


H = 0.005 


TF 


fcis = 3, fc' 16 = 1, fc'/g = 2.9, Ji 5 = 0.01, Ji 6 = 0.01 


kwee 


Cee = 0-115, &^, ee = 1, Viwee = 1-45, V aW ee = 0.25, 
Jawee ' 0.01, Jiwee 0.01 


k 2 5 


fc 25 = 0.01, fc£ 6 = 1, Fi25 = 0.25, K^ee = 0.36, J a25 = 0.01, 

J iwee = 0.01, J i25 = 0.01 


MPF 


k'{ 7 = 0.69, fci 7 = 1.5, fc' 17 = 1.3, k'{ 7 = 1.5, fc'"17 = 1.5 


Trimcr 


fcis = 0.441, k' 18 = 0.882 


a 


k' 19 = 1.5, fc'/g = 0.147, tf diflS = 0.001 



TABLE I: Parameter values for the rescaled system of differential equations 
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TABLE II: Boolean rules for variable CdciZ. 
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FIG. 2: Interaction network of the Boolean model 



Green links correspond to Tik=+1 and red links to ones to Tik=-1 
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TABLE III: Temporal evolution of protein states for the cell-cycle control network. 
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FIG. 3: a): Dynamics of IEP (red curve) and Ste9 (black curve) in the numerical simulation of the system of 
differential equations. IEP (red curve) and 5**69 (black curve), b) Numerical simulation of stationary states for IEP 
(red curve) and Ste9 (black curve), c) Boolean networks of dynamical sequence for IEP (red curve) and Ste9 (black 
curve) . 



